In [2]:
from opnexrad import opnexrad
from opaod import *
from opaodbase import opaod_base
from opecmwf import *
from oppopdens import *
from oplandcover import *
from opsoil import *
from opglim import *
from oplith4 import *
from opgebco import *
from opsa import *
from opmodel import *

from mpl_toolkits.axes_grid1 import make_axes_locatable
import xarray as xr
import os
from cartopy.feature import ShapelyFeature
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import xarray as xr
import cartopy.crs as ccrs
import pyart
import warnings
warnings.simplefilter('ignore')
import time
import importlib
import joblib
%load_ext autoreload
%autoreload 2

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
In [3]:
base = "/home/xiaohe/Documents"
fname = os.path.join(base,'NEXRAD/AirNow/US_Boundary/USBoundary.shp')

shape_feature = ShapelyFeature(Reader(fname).geometries(),
                                ccrs.PlateCarree(), facecolor='none')
In [4]:
AOD_UTC = "2020-01-06T18:30"


mdl_type ="nonex"     #"nex" , "nonex" both nexrad and ecmwf, 

if mdl_type == "nonex":
    print("Load Mdl")
    Mdl = joblib.load(os.path.join(base,"NEXRAD/Model/ResidualAnalysis/US_AODDQF_SA_Z01_ET_Bayes_testMdl/ETreg_ECMWF_random.joblib"), mmap_mode=None)
    predictors = [
           'uv10', 'd2m', 't2m', 'lai_hv', 'lai_lv', 'sp', 'sro', 'tp', 
           'popden', 'landcover', 'soil', 'glim', 'lith4', 'gebco','AOD','DQF','LZA','SZA','SAA']
    
    ### AOD base  ################################################################################
    ds_GOES = xr.open_dataset("noaa-goes16/ABI-L2-AODC/2020/006/18/OR_ABI-L2-AODC-M6_G16_s20200061831177_e20200061833550_c20200061840213.nc",decode_times=False)
    print("trim and grid AOD, save to df")
    df_NG = opaod_base(ds_GOES,plotvar="yes",shape_feature=shape_feature)
    full_extent = df_NG.iloc[:,0:2] # used as base for predicted model merge
    filter_NG = df_NG.dropna()
    
    
    
elif mdl_type == "nex":
    print("Load Mdl")
    Mdl = joblib.load(os.path.join(base,"NEXRAD/Model/ResidualAnalysis/US_AODDQF_SA_Z01_ET_Bayes_testMdl/ETreg_random.joblib"), mmap_mode=None)


    predictors = [
           'z01reflectivity', 'z01velocity',
           'z01spectrum_width', 'z01differential_phase',
           'z01differential_reflectivity', 'z01cross_correlation_ratio', 
           'uv10', 'd2m', 't2m', 'lai_hv', 'lai_lv', 'sp', 'sro', 'tp', 
           'popden', 'landcover', 'soil', 'glim', 'lith4', 'gebco','AOD','DQF','LZA','SZA','SAA']

    ## Read NEXRAD file and retrived iregular 2d coordinates
    grid = pyart.io.read_grid(os.path.join(base,"NEXRAD/Variables/TestingNEXRAD/2020-01-06_183000.nc"))
    NEXZ01, full_extent = opnexrad(grid,shape_feature,plotvar="reflectivity")
    
    ## AOD match ################################################################################
    stamp = "2020-01-06T18:00"
    filter_N = NEXZ01.dropna() # if nexrad is included


    start_time = time.time()
    nc_path = match_goes(filter_N, stamp)
    filter_NG = filter_N.dropna()
    print("-------%s seconds -------" %(time.time() - start_time))
    print("potting AOD")
    aod_plot(nc_path,shape_feature,"AOD")

Load Mdl
trim and grid AOD, save to df
-------3.1011202335357666 seconds -------
In [9]:
### ECMWF #######################################################################################

ds = xr.load_dataset(os.path.join(base,'NEXRAD/Model/Operation/era5_US_Land_2020_1_6_1800.grib'), engine='cfgrib')
match_ecmwf(filter_NG,ECMWF_xr=ds)
filter_NEG = filter_NG.dropna()

# plot emcwf
plot_2decmwf(ds,shape_feature)

retrieving u10 values from ECMWF grid file ......
There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF
retrieving v10 values from ECMWF grid file ......
There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF
retrieving d2m values from ECMWF grid file ......
There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF
retrieving t2m values from ECMWF grid file ......
There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF
retrieving lai_hv values from ECMWF grid file ......
There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF
retrieving lai_lv values from ECMWF grid file ......
There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF
retrieving sp values from ECMWF grid file ......
There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF
retrieving sro values from ECMWF grid file ......
There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF
retrieving tp values from ECMWF grid file ......
There are 20849 grid macthed with the ECMWF, 22212 grid are out of coverage from ECMWF
In [6]:
### popdens  ####################################################################################

start_time = time.time()
usda = combinetiff("PopDens/us_pop15.tif",
                  "PopDens/us_pop20.tif")
# only support pop from 15 to 23
match_popdens(filter_NEG,POPDEN_xr=usda, year="2020")
print("-------%s seconds -------" %(time.time() - start_time))

# ## plot pop dens tif
da20 = xr.open_rasterio("PopDens/us_pop20.tif")
fig = plt.figure(figsize=(20,13))
ax = fig.add_subplot(111)
im =ax.imshow(da20.variable.data[0],vmin=0,vmax=100)
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "5%", pad="3%")
plt.colorbar(im, cax=cax)
plt.tight_layout()
plt.savefig("Popden2020.png")
filter_NEG_P = filter_NEG.dropna()

retrieving population density values in 2020 from CEDAC geotiff file ......
There are 20849 sites macthed with the population density raster, 0 sites are out of coverage from Radar
-------36.47853183746338 seconds -------
In [10]:
#landcover##########################################################################################

start_time = time.time()
da = xr.open_rasterio("LandCover/Reclass_NLCD_Re1kmMajorit.tif")
match_landcover(filter_NEG_P,da)
print("-------%s seconds -------" %(time.time() - start_time))

## plot
lc = xr.open_rasterio("LandCover/Reclass_NLCD_Re1kmMajorit.tif")
fig = plt.figure(figsize=(20,13))
ax = fig.add_subplot(111)
im =ax.imshow(lc.variable.data[0],vmin=1,vmax=18)
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "5%", pad="3%")
plt.colorbar(im, cax=cax)
plt.tight_layout()
plt.savefig("LandCover.png")

filter_NEG_PL = filter_NEG_P.dropna()



### soil ##########################################################################################

start_time = time.time()
da = xr.open_rasterio("Soil/so2015v2.tif")
print("Triming to US extent")
da = da.where(da.y<50,drop=True).where(da.y>25,drop=True).where(da.x>-125,drop=True).where(da.x<-65,drop=True)
match_soil(filter_NEG_PL,da)
print("-------%s seconds -------" %(time.time() - start_time))

filter_NEG_PLS = filter_NEG_PL.dropna()

## plot
fig = plt.figure(figsize=(20,13))
ax = fig.add_subplot(111)
im =ax.imshow(da.variable.data[0],vmin=1,vmax=100)
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "5%", pad="3%")
plt.colorbar(im, cax=cax)
plt.tight_layout()
plt.savefig("soil.png")



### glim ##########################################################################################

start_time = time.time()
da = xr.open_rasterio("GLiM/glim_wgs84_0point5deg.tif")
print("Triming to US extent")
da = da.where(da.y<50,drop=True).where(da.y>25,drop=True).where(da.x>-125,drop=True).where(da.x<-65,drop=True)
match_glim(filter_NEG_PLS,da)
print("-------%s seconds -------" %(time.time() - start_time))

filter_NEG_PLSG = filter_NEG_PLS.dropna()

## plot
fig = plt.figure(figsize=(20,13))
ax = fig.add_subplot(111)
im =ax.imshow(da.variable.data[0],vmin=1,vmax=13)
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "5%", pad="3%")
plt.colorbar(im, cax=cax)
plt.tight_layout()
plt.savefig("glim.png")



### lith4 ##########################################################################################

star_time = time.time()
da = xr.open_rasterio("Lithology4Class/fullareagmnarocktypedecde.tif")
match_lith4(filter_NEG_PLSG, da)
print("-------%s seconds -------" %(time.time() - start_time))

filter_NEG_PLSGL = filter_NEG_PLSG.dropna()

## plot
fig = plt.figure(figsize=(20,13))
ax = fig.add_subplot(111)
im =ax.imshow(da.variable.data[0],vmin=1,vmax=4)
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "5%", pad="3%")
plt.colorbar(im, cax=cax)
plt.tight_layout()
plt.savefig("lith4.png")



### gebco ##########################################################################################

start_time = time.time()
da = xr.open_rasterio("GEBCO_2020_tif/gebco_2020_n50.0_s25.0_w-125.0_e-65.0.tif")
match_gebco(filter_NEG_PLSGL, da)
print("-------%s seconds -------" %(time.time() - start_time))

filter_NEG_PLSGLG = filter_NEG_PLSGL.dropna()



## plot
fig = plt.figure(figsize=(20,13))
ax = fig.add_subplot(111)
im =ax.imshow(da.variable.data[0],vmin=-100,vmax=3300)
divider = make_axes_locatable(plt.gca())
cax = divider.append_axes("right", "5%", pad="3%")
plt.colorbar(im, cax=cax)
plt.tight_layout()
plt.savefig("gebco.png")

retrieving landcover values from  geotiff file ......
There are 20849 sites macthed with the landcover raster, 0 sites are out of coverage from Radar
-------48.6298770904541 seconds -------
Triming to US extent
retrieving soild values from  geotiff file ......
There are 20849 sites macthed with the landcover raster, 0 sites are out of coverage from Radar
-------40.850067138671875 seconds -------
Triming to US extent
retrieving soild values from  geotiff file ......
There are 20849 sites macthed with the landcover raster, 0 sites are out of coverage from Radar
-------40.30678725242615 seconds -------
retrieving values from  geotiff file ......
There are 20849 sites matched with the raster, 0 sites are out of coverage
-------90.83371686935425 seconds -------
retrieving values from  geotiff file ......
There are 20849 sites matched with the raster, 0 sites are out of coverage
-------54.32755208015442 seconds -------
In [11]:
### SA ##########################################################################################

match_sa(filter_NEG_PLSGLG,UTC = AOD_UTC)
filter_NEG_PLSGLGS = filter_NEG_PLSGLG.dropna()
# plot_sa(filter_NEG_PLSGLGS,shape_feature,plotvar="SZA")


### Model ##########################################################################################
modelPredict(filter_NEG_PLSGLGS, full_extent, predictors, Mdl,shape_feature)

Calculating LZA
Calculating SA 
Progress: [----------------------------->] 100 %
Out[11]:
Latitude Longitude predictions
0 25.0 -125.000000 NaN
1 25.0 -124.874739 NaN
2 25.0 -124.749478 NaN
3 25.0 -124.624217 NaN
4 25.0 -124.498956 NaN
... ... ... ...
158395 50.0 -65.501044 NaN
158396 50.0 -65.375783 NaN
158397 50.0 -65.250522 NaN
158398 50.0 -65.125261 NaN
158399 50.0 -65.000000 NaN

158400 rows × 3 columns

In [ ]:
In [ ]:
In [ ]:
In [ ]:
# # use 2km resoluation as base grid, need to regrid popdens landcover lith4 gebco to 2km grid

# da20 = xr.open_rasterio("PopDens/us_pop20.tif")
# ds = da20

# xx, yy = np.meshgrid(ds.x.data, ds.y.data)

# x_flat = xx.flatten()
# y_flat = yy.flatten()
# x_flat[x_flat<-125] = None
# x_flat[x_flat>-65] = None
# y_flat[y_flat>50] = None
# y_flat[y_flat<25] = None


# arr = ds.data
# arr_flat = arr.flatten()

# trimmed = pd.DataFrame()
# trimmed["Latitude"] = y_flat
# trimmed["Longitude"] = x_flat
# trimmed["popdens"] = arr_flat

# ## filter out None values which are outside US
# trimmed_drop = trimmed.dropna(subset=["Latitude"])
# trimmed_drop = trimmed_drop.dropna(subset=["Longitude"])

# ## 10km grid which is the same as nexrad extent
# xi = np.linspace(-125,-65,2400)
# yi = np.linspace(25,50,1650)
# xi,yi = np.meshgrid(xi,yi)

# zi = griddata((trimmed_drop["Longitude"], trimmed_drop["Latitude"]),trimmed_drop["popdens"], (xi,yi), method='linear')

# ds_NEXREG = xr.DataArray(zi, coords=[yi[:,1], xi[0]], dims=["lat", "lon"])

# fig = plt.figure(figsize=[20,13])
# ax = plt.subplot(projection=ccrs.PlateCarree())
# ax.add_feature(shape_feature,edgecolor='blue')
# ax.gridlines(draw_labels=True)
# ax.coastlines()
# ds_NEXREG.plot(cmap=plt.cm.Reds, vmin=0,vmax=100, x='lon', y='lat')
# plt.savefig("poptestxrlinear.png")